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A  HIGH-RESOLUTION  TARGET-TRACKING  CONCEPT 
USING  SPECTRAL  ESTIMATION  TECHNIQUES 


1.  INTRODUCTION 

The  tracking  of  targets  in  the  presence  of  nearby  strong  interference  sources  is  a  problem  area  of 
considerable  interest  in  radar  systems,  primarily  because  conventional  tracking  radars  may  experience 
large  errors,  track-breaking,  or  complete  disruption.  Some  early  proposed  solutions  to  this  problem 
evolved  frrtr.  the  growing  adaptive  array  antenna  technology  of  the  1970’s.  For  example,  a  paper  by 
White  [1]  discusses  the  hish-resolution  aspect  in  considerable  detail.  He  presents  the  problem  of  track¬ 
ing  radar  targets  in  the  low-angle  regime  where  conventional  tracking  radars  encounter  difficulty 
because  of  the  presence  of  a  strong  surface-reflected  ray.  Starting  with  a  classical  maximum-likelihood 
analysis  of  the  problem  of  two  closely  spaced  targets,  he  develops  two  techniques  which  are  theoreti¬ 
cally  capable  of  dealing  with  the  multipath  problem.  A  more  recent  contribution  in  the  multipath  area 
is  the  work  of  Cantrell  et  at.  [2],  who  proposed  a  simple,  closed  form  solution  for  an  array  divided  into 
three  subapertures.  The  accuracy  of  their  technique  compares  favorably  with  the  maximum  likelihood 
estimate  which  uses  all  of  the  individual  elements. 

The  first  extension  of  fully  adaptive  arrays  to  angle  estimation  in  external  noise  fields  is  the  con¬ 
tribution  of  Davis  et  al.  [21,  who  developed  an  algorithm  based  on  the  outputs  of  adaptively  distorted 
sum  and  difference  beams.  The  adaptive  beams  filter  (null)  the  external  noise  sources,  and  distortion 
correction  is  then  applied  in  the  resultant  monopulse  output  angle  estimate  for  a  selected  tracking  beam 
direction. 

In  recent  years,  a  closely  related  new  technology  which  derives  largely  from  modern  spectral  esti¬ 
mation  theory  has  been  emerging,  and  it  appears  to  offer  additional  opportunities  for  achieving 
significant  improvements  in  the  detection/tracking  of  closely  spaced  multiple  sources  and  targets. 
In*erest  is  stimulated  because  of  reported  resolution  capabilities  which  extend  well  beyond  the  conven¬ 
tional  windowed  Fourier  transform.  Gabriel  (4,51  has  emphasized  two  processing-domain  performance 
characteristics  for  these  techniques  which  apply  to  the  tracking  problem: 

•  "superresolution"  which  implies  the  resolution  of  two  or  more  source.®  within  the  conven¬ 
tional  beamwidth  of  an  RF  antenna  aperture,  and 

•  "absence  of  sidelobes"  which  implies  the  resolution  of  two  or  more  sources  of  unequal 
strengths  when  those  sources  are  spaced  more  than  a  beamwidth  apart. 

(Quotation  marks  are  used  to  avoid  misinterpretation  of  the  conventional  radar  system  definitions  of 
resolution  and  sidelobes.)  These  characteristics  are  demonstrated  via  simulation  examples  for  several 
spectral  estimation  algorithms  including  linear  prediction  methods,  the  maximum  likelihood  method, 
and  eigenvalue/eigenvector  analysis  techniques.  References  4  and  5  also  address  the  difficult  radar 
coherent-source  problem  area  which  involves  a  phase-dependent  signal-to-noise  ratio  (SNR)  penalty. 

A  contribution  which  proposes  both  adaptive  processing  and  high-resoiution  spectral  estimation  is 
that  of  Chapman  et  al.  (6),  who  discusses  the  mainbeam  jamming  problem.  Two  techniques  were 
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investigated.  First,  a  linear  prediction  extrapolation  procedure  was  employed  to  extend  the  apparent 
aperture  of  an  antenna  array  and  thereby  achieve  a  narrow  notch  on  the  jammer  while  providing  gain  in 
the  look  direction  of  the  target.  Second,  the  maximum  likelihood  procedure  (Section  4.4)  was  used  to 
provide  adaptive  filtering  and  estimate  target  angle  via  beam  scan. 

The  remainder  of  this  report  is  divided  into  three  principal  sections.  Section  2  discusses  the 
overall  tracking  system  concept  and  references  the  sections  containing  technical  details.  Section  3  con¬ 
siders  several  tracking  performance  examples  which  are  based. on  computer  simulation  snapshot  data, 
and  Section  4  discusses  the  details  of  the  associated  algorithm/processing  theory,  analysis,  and  develop¬ 
ment 

2.  TRACKING  SYSTEM  CONCEPT 

An  all-digital  processing  system  is  defined  which  commences  with  complex  data  sampling  of  the 
array  element  signals  (or  beam-port  signals  if  it  is  a  beamformer  type  of  antenna).  It  is  assumed  that 
conventional  receiver  techniques  would  be  employed  prior  to  the  signal  analogue-to-digital  (A/D)  con¬ 
verters  as  illustrated  in  Fig.  1;  i.e.,  the  system  may  include  transmit-receive  (TR)  devices,  low-noise  RF 
amplifiers,  local  oscillator  (LO)  mixers,  intermediate  frequency  (IF)  amplifiers,  and  synchronous 
baseband  detection  with  in-phase  and  quadrature  (I&Q)  video  output  to  sample  and  hold  (S/H)  circuits 
followed  by  the  A/D  converters.  Snapshots  would  nominally  occur  at  the  Nyquist  sampling  rate 
corresponding  to  the  video  bandwidth,  so  that  a  radar-oriented  person  may  view  them  as  range  bin  time 
samplings.  For  this  interim  report,  actual  measured  snapshot  data  were  not  used  in  tracking  perfor¬ 
mance  evaluations.  Rather,  a  snapshot  signal  model  described  in  Section  4.1  was  used  to  generate  the 
complex  digital  signal  data  via  computer  simulation. 


Fig.  1  —  Typical  RF  receiver  techniques  associated 
with  A/D  complex-data  sampling 


Figure  2  shows  the  system  concept  for  processing  the  digital  signals.  Starting  on  the  left-hand 
side,  the  system  continuously  computes/updates  a  sample  covariance  matrix  R.  A  number  of  tech¬ 
niques  are  available  for  computing  R,  and  four  of  these  are  described  in  Section  4.2.  Of  particular 
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Fig.  2  —  Adaptive  array  tracking  system  concept 


significance  is  the  fact  that  ft  may  to  dimensioned  either  equal  to  or  less  than  the  number  of  array  ele¬ 
ments;  i.e.,  the  model  order  of  the  estimate  is  selectable  per  subaperture  averaging  option  choice. 

OfF-iine  processing  on  ft  is  then  conducted  at  periodic  intervals  to  estimate  the  locations  and  rela¬ 
tive  power  levels  of  interference  sources  via  a  catalogue  of  spectral  estimation  algorithms.  Literally 
dozens  of  algorithms  are  now  available  in  the  literature,  ranging  from  fast  simple  routines  [4,7]  for 
estimating  strong  sources,  to  the  more  sophisticated  eigenanaiysis  routines  for  estimating  weak  sources. 
Because  of  the  growing  interest  and  superior  performance  of  eigenvalue/eigenvector  decomposition 
techniques,  Sections  4.3  to  4.5  are  devoted  to  a  detailed  discussion/analysis  with  pertinent  references. 
Desirable  features  of  any  algorithm  are  robustness,  fast  root-finding,  and  the  capability  of  solving  for 
relative  source  strength.  This  latter  feature  is  important  in  source  classification  and  crucial  for  screen¬ 
ing  out  "false  alarms." 

Having  estimated  the  locations  and  relative  power  levels  of  interference  sources,  the  central  pro¬ 
cessor  unit  (CPU)  then  applies  these  data  to  the  computation  of  optimized  adaptive  spatial  filter 
weights.  These  may  be  developed  via  a  number  of  different  schemes,  and  Section  4.7  discusses  two 
approaches  with  an  option  for  "softened"  spatial  filtering.  The  question  may  be  raised  as  to  why  the 
filter  weights  are  not  computed/developed  directly  from  the  data  samples  as  in  a  conventional  adaptive 
array.  A  number  of  reasons  exist  for  separating  source  estimation  from  filter  weight  computation. 

1.  Source  estimation  often  requires  flexible  subaperture  techniques,  i.e.,  a  model  order 
which  is  less  than  full  aperture.  (See  Section  4.2.4.) 
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2.  Search/track  detection  usually  demands  use  of  the  full  aperture  at  all- times  because  of 
SNR  and  resolution  considerations. 

3.  It  is  easier  to  estimate  and  implement  prewhitening  filtering  to  remove  slowly  changing 
interference  sources  and  colored  noise  distributions. 

4.  It  is  easier  to  modify  the  filtering  for  purposes  other  than  source  cancellation,  i.e.,  to 
reduce  distortion  effects  cr  to  introduce  time-modulated  weights. 

5.  One  can  readily  add  real-time  adaptive  filter  weight  perturbations  to  permit  fast  response. 

6.  Certain  constrained,  partially  adaptive  systems  are  unable  to  assign  their  degrees  of  free¬ 
dom  automatically,  for  example,  a  beamspace  adaptive  antenna  system  with  a  low  sidelobe 
constraint  and  only  a  few  beams  available. 

The  right-hand  side  of  Fig.  2  shows  that  the  CPU  also  has  control  over  a  fast-memory  storage 
capability  which  could  be  used  in  several  different  ways  including 

•  dump-storage  of  a  block  of  contiguous  data  snapshots;  for  example,  one  may  store  the 
snapshots  associated  with  a  particular  sample  covariance  matrix,  for  later  retrieval  when 
that  filter  update  occurs; 

•  time-gated  snapshot  storage  based  on  interference  source  classification,  for  later  selective 
retrieval; 

•  time-gated  snapshot  storage  based  on  target-search  considerations;  and 

•  "bucket-brigade  bumping”  snapshot  storage  to  obtain  a  selected  time  delay  of  the 
snapshots. 

Thus  one  may  select  either  the  current  snapshot  or  a  storage  snapshot  for  feeding  into  the  filter  weights 
which  the  CPU  updates. 

Next,  the  filtered  signal  output  residue  is  fed  into  a  beamformer  which  is  weighted  to  produce  the 
desired  search  and  monopulse  track  beams  for  target  detection/tracking.  Section  4.8  discusses  the  un¬ 
avoidable  problem  of  target  wavefront  distortion  which  can  cause  serious  tracking  errors  if  not  dealt 
with  and  points  out  two  distortion  correction  procedures. 

Recognize  that  Fig.  2  is  intended  to  represent  a  concept  rather  than  a  specific  system  design.  Any 
particular  system  design  application  would  undoubtedly  require  variants  to  achieve  the  most  efficient 
processor  implementation.  The  concept  is  summarized  in  the  following  five  steps: 

1.  Convert  A/D  at  elements/beam  ports. 

2.  Compute/update  sample  covariance  matrix  R. 

•  model  order  selectable  (need  not  be  known) 

•  subaperture  technique  incorporated 

3.  Process  R  off-line  periodically 

•  by  using  fast  inverse  routines  to  estimate  strong  sources, 

4 


.  W.  r  . 


W.V.V  V.V.v.V 


.  V.  • 


NRL  REPORT  8797 


•  by  using  eigenanaly.  is  routines  to  estimate  weak  sources, 

•  by  rootfinding  on  optimum  weights  to  pick  out  source  locations; 

•  then  compute  source  power  matrix; 

•  theu  screen  false  alarms. 

4.  Compute  filter  weights  for  full  aperture  from  updated  source  information. 

5.  Apply  conventional  tracking  to  the  filtered  snapshots,  incorporating  wavefront  distortion 
correction. 

3.  TRACK  PERFORMANCE  EXAMPLES 

To  demonstrate  the  concept,  a  series  of  simple  radar-simulation  examples  are  given  in  this  section 
wherein  a  single  target  appears  in  a  single  range  bin  per  pulse  repetition  frequency  (PRF)  period,  and 
one  to  several  interference  sources  appear  in  N  adjacent  range  bins  per  PRF  period.  All  snapshots  were 
simulated  from  the  model  described  in  Section  4.1,  by  using  an  8-element  linear  array  with  half¬ 
wavelength  spacing  as  the  receive  aperture.  The  examples  are  identified  as  follows: 

Case  C8L6JVG0B  -  six  sources, 

Case  C8L1JXX  -  one  source,  | 

Case  C8L2J3C  -  two  sources. 

Case  C8L3J3C  -  three  sources,  and 

Case  C8L3J7B  •  three  sources  plus  spatial  noise. 

3.1  Case  C8L6JVC*B  -  Six  Sources 

This  case  was  originally  designed  in  January  1981  for  an  8-element  array  for  the  purpose  of  stress¬ 
ing  spectral  estimation  algorithms.  Figure  3  illustrates  its  range-azimuth  geometry,  which  consists  of  a 
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single  3  dB  SNR  target  located  at  -4°  azimuth  and  appearing  in  a  single  range  bin  per  PRF  period.  Six 
interference  sources  appear  in  N  adjacent  range  bins  per  PRF  period.  Two  of  these  sources  are  10  dB 
SNR,  95%  correlated,  located  at  —54°  and  —42°  (0.57  beamwidth  apart);  and  four  are  noncoherent 
sources  of  —9,  30,  10,  and  20  dB  SNR,  located  at  0s,  15°,  22°,  and  30°  respectively.  The  correlated 
sources,  the  close  spacings,  the  large  dynamic  range  of  source  strengths,  and  the  consumption  of  six 
out  of  seven  degrees  of  freedom  causes  difficulty  for  most  spectral  estimation  algorithms,  as  shown  in 
Fig.  4.  Note  in  Fig.  4(a)  that  the  conventional  beam  scan  (windowed  Fourier  transform)  is  hopeless,  as 
expected,  because  of  its  beamwidth  and  sidelobe  limitations.  But  also  note  that  the  SLC  (sidelobe  can¬ 
celler  algorithm  similar  to  least-squares  linear  prediction)  has  failed  to  resolve  the  highly  correlated 
sources  and  the  two  weaker  noncoherent  sources  located  at  0®  and  22°.  In  Fig.  4(b)  the  maximum 
Likelihood  method  (MLM)  algorithm  also  performs  poorly,  but  the  multiple  signal  classification 
(MUSIC)  algorithm  (Section  4.4.2)  was  able  to  locate  all  six  sources  correctly  via  its  noise  eigenvectors. 
The  point  to  be  made  is  that  this  is  not  an  easy  case  to  resolve. 

This  first  case  will  be  presented  in  greater  walk-through  detail  than  the  others  to  familiarize  the 
reader  with  the  sequence.  Starting  with  the  computation  of  R,  the  simple  block  averaging  method  (Sec¬ 
tion  4.2.1)  was  used  with  N  ™  256  snapshots.  Further,  the  filter  weight  update  period  was  selected  to 
coincide,  so  that  a  new  copy  of  R  was  fetched  for  processing  at  the  end  of  each  256  snapshots  period. 
Two  spectral  estimation  algorithms  were  applied  to  R:  the  SLC  algorithm  and  the  principal  eigenvector 
Gram-Schmidt  (PEGS)  algorithm  (Section  4.4.3).  The  SLC  algorithm  develops  the  Wiener  optimum 
weight  [4,5],  per  Eq.  (25)  in  Section  4.4.1,  W„  -  /aR_1S®,  where  S®'«*  [0,  0,  0,  ....  1].  Rootfinding 
(Section  4.6)  on  the  first  four  updates  of  W„  resulu  in  the  Z-piane  roots  plots  shown  in  fig.  5(a).  Note 
that  only  three  roots  are  close  to  the  unit  circle  and  would  be  culled,  leading  then  to  the  corresponding 
digital  power/location  estimates  shown  in  Fig.  5(b).  Performance  comments  here  are  the  same  as  those 
associated  with  Fig.  4(a);  i.e.,  the  SLC  algorithm  has  failed  to  resolve  the  highly  correlated  sources  and 
the  two  weaker  noncoherent  sources  located  at  0®  and  22®. 

The  PEGS  eigenanalysis  algorithm  applied  to  the  same  four  R  update  trials  results  in  the  roots 
and  digital  power/location  estimates  shown  in  Fig.  6.  Note  that  the  estimates  are  now  reasonably  close 
to  the  truth  given  in  Fig.  3,  with  the  greatest  variance  exhibited  by  the  highly  correlated  sources  and 
the  weak  —9  dB  source. 

After  the  source  data  are  estimated,  the  adaptive  spatial  filter  weights  are  updated  accordingly 
(Section  4.7),  and  we  can  then  commute  the  associated  spatial  filter  insertion  loss.  Figure  7  illustrates 
the  results  for  the  four  trials  in  this  example,  computed  from  Eq.  (50)  for  a  —  1.  Insertion  loss  plots 
are  very  helpful  for  visualizing  certain  target/interference-source  relationships. 

•  Detection/tracking  is  most  favorable  in  those  regions  where  the  insertion  loss  is  low. 

•  Insertion  loss  cost  against  a  desired  target  increases  rapidly  as  the  target  gets  close  to  a 
strong  interference  source. 

•  If  the  target  strength  is  sufficient  to  overcome  the  insertion  loss,  then  tracking  within  a 
beamwidth  of  a  strong  interference  source  should  be  feasible. 

•  Insertion  loss  and  tracking  beam  distortion  are  related. 
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(a)  Conventional/SLC  algorithm 


(b)  MLM/eigenvector  algorithm,  aymmetric  matrix  constraint 

Fig.  4  —  Spatial  spectrum  estimate  for  six  sources:  two  10  dB  sources  95%  correlated  at 
-54*  and  -42*  (0.57  beamwidth  apart);  four  noncoherent  sources  of  -9,  30,  10,  and 
20  dB  strengths  located  at  0*,  15*,  22*,  and  30*;  8-element  linear  array,  sample  covari¬ 
ance  matrix  of  1024  snapshots 


U 

H 

/ 

V 


V 
■J* 

V 


'A 


S 

!*, 


S 


f 


*> 


S 

•% 

A 


S 


i 


.  ,  •  a  .  v'  --  /■  .*  .*■  *VI  'A  *  •  **>  *V  V  ’.*  *.'■  *, 


7 


WILLIAM  F.  GABRIEL 


SPATIAL  AMQLE  IM  DEGREES 
(b)  Digital  relative  power  vi  location  anile 


Fif.  6  —  Spatial  apectnim  estimates  Tor  Case  CiL6JVC^B,  derived  from  PEGS  eitenanatyiis 
algorithm  optimum  weights,  four  trials  of  256  snapshots 
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SPATIAL  ANGLE  IN  DEGREES 
Fig.  7  —  Adaptive  spatial  filter  insertion  loss 


Section  4.8  discusses  the  formation  of  the  tracking  beams  and  their  resultant  equivalent  distortion 
due  to  the  spatial  filtering.  Figure  8  summarizes  the  particular  distortions  associated  with  this  case; 
these  distortions  are  computed  for  tracking  beams  with  a  boresight  pointing  angle  of  0°.  The  reader  will 
benefit  from  comparing  these  patterns  against  the  sidelobe  source  case  of  Fig.  10.  Figure  8(a)  shows 
the  distorted  adapted  left  and  right  tracking  beams,  where  the  right  beam  is  particularly  adversely 
affected  because  of  its  overlap  with  the  three  strong  sources  located  at  15,  22,  and  30°.  Figure  8(b) 
shows  the  resultant,  distorted,  adapted  sum  and  difference  tracking  beams  (a  plot  of  W,  and  \\d  as 
given  in  Eqs.  (65,  66)),  where  we  note  that  the  sum  beam  has  been  shoved  over  to  the  left  and  the 
difference  beam  has  almost  lost  its  right  positive  lobe.  Figure  8(c)  is  the  associated  spatial  filter  inser¬ 
tion  loss  over  the  tracking  beam  region,  an  expanded  portion  of  Fig.  7.  And,  Fig.  8(d)  illustrates  the 
associated  track  angle  estimate  (A/X)  distortion,  computed  from  Eqs.  (64,  65,  66).  Note  the  cusp 
which  occurs  at  the  location  of  the  strong  30  dB  source.  This  is  undesirable,  of  course,  since  it  could 
result  in  ambiguous  track  angle  estimates  in  that  vicinity.  However,  the  insertion  loss  plot  tells  us  that 
it  would  be  virtually  impossible  to  see  any  target  return  from  that  vicinity  and,  therefore,  one  can  disre¬ 
gard  it. 

A  verification  data  point  is  included  for  our  true  target  location  of  —4*.  It  was  obtained  from 
monopulse  track  angle  estimates  averaged  over  16  trials.  Figure  9  shows  a  representative  plot  of  track¬ 
ing  angle  estimate  vs  time  in  PRF  periods,  where  an  averaging  factor  of  10  is  incorporated  into  the  plot 
outputs  to  reduce  the  dispersion  of  the  plot  and  to  keep  the  various  curves  separated  reasonably  well. 
For  example,  each  track  angle  estimate  plot  point  is  an  average  of  the  previous  10  monopulse  track 
angle  estimates  computed  for  the  range-gated  target  pulses.  The  +  and  —  symbols  represent  X  and  A 
beam  output  samplings  just  prior  to  the  target  pulse.  These  samplings  verify  that  the  interference 
sources  have  been  attenuated  down  to  receiver  noise  level  at  the  output;  i.e.,  they  should  merge  into 
random  receiver  noise  with  a  mean  of  0  dB. 
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Z.  TARGET  SfeR  ■  11  dB 
INTEGRATION  -  10  PULSES  NONCOHERENT 
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Fig.  9  —  Time-tequence  plots  of  sum  besm  output,  difference  beam  output,  and 
mo  no  pulse  tracking  angle  estimate  output  for  range-gated  target  pulses,  10-pulse  non¬ 
coherent  integration 


3.2  Case  C8L1JXX  -  One  Source  . 

/ 

These  cases  were  designed  to  demonstrate  how  the  track  estimate,  distortion  js  affected  by  the 
position  of,  the  tracking  beams  with  respect  to  the  interference  soured  A  single  7  dB  SNR  target 
appears  in  a  single  range  bin  per  PRF  period  and  is  located  within  ±  1  beamwidth  of  boresight  Table  1 
contains  the  case  identification  numbers  corresponding  to  the  locations  of  a  single  10  dB  SNR  interfer¬ 
ence  source  which  appears  in  all  range  bins  per  PRF  period.  The  tracking  beams  of  our  8  element  array 
have  a  boresight  pointing  angle  of  0#,  and  a  uniform  illumination  beamwidth  B  —  14.5".  \ 


Table  1  —  Case  Identification 
Numbers 


Case 

Number 


C8L1J8X 

C8L1J6X 

C8L1J4X 

C8L1J3X 

C8L1J2X 

C8L1J0X 


Source  Location 

9 

(9jB) 

-30° 

-2.0 

-22° 

-1.5 

-14.5* 

-1.0 

-10.81° 

-0.75 

-  7.2° 

-0.50 

0° 

0 

The  simple  block  averaging  method  was  used  to  compute  R,  with  N  —  480  snapshots,  and  the 
PEGS  eigenanalysis  algorithm  was  applied.  Estimating  a  single  source  is  almost  trivial,  and,  as  a  result, 
showing  the  roots  or  the  digital  power/location  estimates  is  not  necessary. 
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One  of  the  general  characteristics  demonstrated  by  this  series  of  simulation  tracking  runs  was 
that  tracking  beam  distortion  was  rather  negligible  so  long  as  the  interference  source  remained  in  the 
sidelobe  regions  of  the  tracking  beams,  which  is  beyond  approximately  ±  22°  for  our  example.  Figure 
10  illustrates  this  point  for  a  source  location  of -2.0  beamwidths.  Thus,  distortion  builds  up  only  when 
the  tracking  main  beams  begin  to  overlap  a  strong  source.  Figure  1 1  illustrates  a  typical  case  where  the 
interference  source  is  located  0.75  beamwidth  from  boresight.  Note  the  severe  distortion  of  the  left 
beam  in  Fig.  11(a)  and  the  resultant  distortion  of  the  sum  and  difference  beams  in  Fig.  11(b).  A 
severe  cusp  occurs  in  the  track  estimate  distortion  of  Fig.  11(d)  at  the  position  of  the  interference 
source,  and  this  could  lead  to  ambiguous  track  angle  estimates  in  that  vicinity.  One  solution  to  this 
ambiguity  is  to  simply  increase  the  filter  insertion  loss  (Fig.  11(c))  where  needed  (similar  to  adding  a 
synthetic  interference  source)  so  as  to  remove  any  possibility  of  a  target  return  getting  through  on  the 
wrong  side.  Figure  11(d)  includes  six  verification  data  points  for  true  target  locations  of  — 1  0,  —0.5, 
-0.25,  0,  +0.5,  and  +1.0  beamwidth.  These  points  were  obtained  from  monopulse  track  angle  esti¬ 
mates  averaged  over  60  trials.  Note  that  satisfactory  tracking  was  obtained  within  0.25  beamwidth  of 
the  interference  source. 

Figure  12  illustrates  maximum  distortion  ambiguity  which  occurs  for  the  special  case  when  track¬ 
ing  beam  boresight  is  directly  steered  onto  a  strong  source.  Tracking  targets  under  this  condition  is  not 
possible,  and  it  obviously  should  be  avoided. 

3.3  Case  C8L2J3C  -  Two  Sources 

The  two-source  cases  were  designed  on  the  basis  of  results  from  the  single  source  cases  to  demon¬ 
strate  how  the  track  estimate  distortion  is  affected  by  the  position  of  the  tracking  beams  with  respect  to 
the  two  sources,  particularly  when  the  tracking  main  beams  overlap  the  sources.  The  example  chosen 
for  presentation  is  like  a  superposition  of  Case  numbers  C8L1J6X  and  C8L1J3X  from  Table  1;  i.e., 
there  are  two  10  dB  SNR  sources  located  at  —1.5  and  -0.75  beamwidths  from  tracking  beam  boresight. 

Again,  the  processing  involved  simple  block  averaging  to  update  R,  with  N  “  480  snapshots;  the 
PEGS  eigenanalysis  algorithm  was  applied,  and  the  estimation  of  the  two  sources  was  easy  because  of 
their  relatively  large  separation  of  0.75  beamwidth. 

Figure  13,  which  illustrates  the  resulting  adaptive  spatial  filter  functions,  may  be  directly  com¬ 
pared  to  Fig.  11  because  of  the  identical  location  of  the  close  source  at  -0.75  beamwidth.  Note  that 
there  is  some  further  degradation  in  all  four  plots  on  the  left-hand  side  because  of  the  presence  of  the 
second  source,  but  that  performance  on  the  useful  right-hand  side  has  not  been  degraded.  In  fact,  the 
track  estimate  distortion  curve  in  Fig.  13(d)  is  now  greatly  improved  because  the  possibility  for  ambi¬ 
guity  in  the  vicinity  of  the  cusp  has  been  removed  by  the  insertion  loss  notch.  Five  verification  data 
points  are  included  for  true  target  locations  of  -0.5,  —0.25,  0,  +0.5,  and  +1.0  beamwidth,  averaged 
over  60  monopulse  estimate  trials.  Again,  satisfactory  tracking  was  obtained  within  0.25  beamwidth  of 
the  close  source,  which  means  that  the  target  plus  both  interference  sources  were  located  within  a 
beamwidth. 

3.4  Case  C8L3J3C  -  Three  Sources 

This  case  was  designed  by  adding  a  third  source  between  the  above  two  sources;  i.e.,  there  are 
now  three  10  d3  SNR  sources  located  at  —1.5,  —1.12,  and  -0.75  beamwidths  from  tracking  beam 
boresight  at  0°.  The  close  spacing  of  0.375  beamwidth  between  the  three  sources  makes  it  a  difficult 
case  to  estimate  accurately;  as  a  result,  more  processing  results  will  be  given  for  this  example.  Simple 
block  averaging  was  again  employed  to  update  R  every  480  snapshots,  and  the  SLC  algorithm  (Eq. 
(25))  and  PEGS  eigenanalysis  algorithm  were  both  applied. 
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SRpTlAl  ANCLE  IN  DECREES  SPATIAL  ANCLE  IN  DEGREES 

00  Adapted  tracking  beams  (c)  Spatial  insertion  loss 


-is  ^  e  is  39  -i  a  t  2 

SPATIAL  ANCLE  IN  DECREES  TRUE  TARCET  ANCLE  IN  IEANUIDTHS 


(b)  Adapted  sum  and  difference  (d)  Treck  estimate  (A/I)  distortion 

Fig.  11  —  Adaptive  processing  spatial  filter  functions  related  to  tracking  angle  estimation  for  Case  C8L1J3A-I 
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(a)  Adapted  tricking  beams 


(c)  Spatial  insertion  loss 


Fig.  13  —  Adaptive  processing  spatial  filter  functions  related  to  tracking  angle  estimation  for  Case  C8L2J3C-I 
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Rootfinding  on  the  first  eight  updates  of  W0  from  the  SLC  algorithm  resulted  in  the  Z-plane  roots 
shown  in  Fig.  14(a).  Note  that  only  two  roots  are  close  to  the  unit  circle,  and  these  roots  correspond  to 
the  two  outer  sources  located  at  -22°  and  -10.8s.  The  SLC  algorithm  completely  fails  to  resolve  the 
third  middle  source. 

By  contrast,  the  PEGS  eigenanalysis  algorithm  applied  to  the  same  eight  R  updates  results  in  the 
roots  shown  plotted  in  Fig.  14(b).  Here,  ali  three  sources  are  readily  resolved,  and  their  roots  hug  the 
unit  circle  closely.  Even  the  off-circle  roots  show  considerably  less  scatter.  The  corresponding  digital 
plot  of  estimated  relative  power  vs  estimated  source  location  is  shown  in  Fig.  15(a),  where  we  note  an 
accurate  estimate  of  the  power  level  of  the  third  middle  source  in  addition  to  its  location.  Figure  15(b) 
then  illustrates  the  resulting  filter  insertion  loss,  computed  for  a  “  1. 

Figure  16  shows  a  summary  of  the  four  adaptive  spatial  filter  functions  for  tracking.  In  compari¬ 
son  to  Figs.  11  and  13,  note  that  only  the  left-hand  side  of  the  four  plots  shows  any  significant  changes, 
and  that  the  right-hand  side  of  Fig.  16(d)  is  well  protected  from  any  cusp  ambiguity  because  of  the 
deep  insertion  loss  notch  extending  over  the  entire  interference  source  region.  Five  verification  data 
points  are  inclrded  for  true  target  locations  of  -0  5,  -0.25,  0,  +0.5,  and  +1.0  beamwidth,  averaged 
over  60  monopulse  estimate  trials.  A  remarkable  feature  here  is  that  satisfactory  tracking  has  been 
obtained  even  with  three  strong  interference  sources  located  within  a  beamwidth  of  the  target,  i.e.,  the 
verification  data  point  for  the  true  target  location  of  -0.5  beamwidth. 

As  an  additional  indicator  of  tracking  results  for  this  example.  Fig.  17  illustrates  a  typical  plot  of 
track  angle  estimate  vs  time  in  PRF  periods.  An  averaging  factor  of  10  is  incorporated  into  the  plot 
outputs  to  reduce  the  dispersion  of  the  plot  points  and  to  keep  the  various  curves  separated  reasonably 
well.  For  example,  each  track  angle  estimate  plot  point  is  an  average  of  the  previous  10  monopulse 
track  angle  estimates  computed  for  the  range-gated  target  pulses.  The  +  and  -  symbols  represent  I 
and  A  beam  output  samplings  just  prior  to  the  target  pulse,  and  they  serve  to  verify  that  the  interfer¬ 
ence  sources  have  been  attenuated  down  to  receiver  noise  level  at  the  output;  Le.,  they  should  merge 
into  random  receiver  noise  with  a  mean  of  0  dB. 

3.5  Case  C8L3J7B  -  Three  Sources  plus  Spatial  Noise 

This  case  is  included  to  illustrate  the  performance  of  the  PEGS  eigenanalysis  algorithm  when  dis¬ 
tributed  background  spatial  noise  (colored  noise)  is  added  to  the  previous  three-source  situation.  The 
spatial  noise  in  the  simulation  consisted  of  35  random  noise  generators  of  —7  dB  SNR  level  distributed 
uniformly  in  sin  0  space  from  —60°  to  0s.  Figure  18  shows  the  spatial  spectrum  estimates  and  the  asso¬ 
ciated  filter  insertion  loss  obtained  for  the  spatial  noise  alone,  for  eight  trials  of  480  snapshots  each. 
Note  that  the  PEGS  algorithm  estimates  the  distributed  noise  sources  quite  well  by  using  five  roots 
(degrees  of  freedom)  to  cover  that  extensive  spatial  region;  i.e.,  it  estimates  on  the  basis  of  its  8- 
element  array-aperture  resolution  capability.  The  insertion  loss  plot  then  provides  the  "envelope"  out¬ 
line  of  the  spatial  block  of  background  noise,  somewhat  similar  to  an  upside-down  version  of  a  win¬ 
dowed  Fourier  transform.  This  information  could  be  used  to  provide  prewhitening  spatial  filtering  of 
the  background  noise. 

The  addition  of  the  three  10  dB  sources  at  -10.8s,  -16.3s,  and  —22s  then  results  in  the  overall 
estimates  shown  in  Fig.  19.  Note  that  the  PEGS  algorithm  continues  to  use  only  five  roots  in  estimat¬ 
ing  the  combined  situation,  but  we  can  readily  discern  the  superposition  of  the  two  sets  of  sources. 
Some  detail  is  lost  since  the  spatial  width  of  the  background  noise  has  shrunk  slightly,  and  the  three  10 
dB  sources  are  estimated  as  two  sources  (Fig.  19(a)).  This  latter  result  occurs  because  the  background 
noise  degrades  the  relative  SNR  of  the  three  sources.  Overall  tracking  distortion  effects  did  not  greatly 
differ  from  Fig.  16. 
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(a)  SLC  algorithm 


(b)  PEGS  ebteianalysis  algorithm 

Fig.  1*  —  Z- plane  roots  derived  from  optimum  weights, 
efebt  trials  of  *80  snapshots  each.  Case  C8L3J3C 
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4.  ALGORITHM  THEORY/ANALYSIS/DEVELOPMENT 

A  considerable  number  of  algorithms  in  the  literature  could  be  applied  to  this  problem  area,  but 
attention  hereiu  is  focused  primarily  on  eigenvalue/eigenvector  decomposition  techniques,  sometimes 
referred  to  as  principal  eigenvector  methods  or  singular  value  decomposition  techniques.  The  reasons 
for  preferring  this  class  of  algorithm  is  that  these  techniques  demonstrate  superior  performance  in  the 
form  of  robustness  against  noise  perturbations,  maximum  detection  sensitivity,  and  maximum  informa¬ 
tion  potential.  Adequate  references  are  given  below  so  that  the  reader  may  gain  a  better  perspective  of 
this  rather  broad  technical  area  and  may  also  consult  more  detailed  presentations  of  the  basic 
theory/analysis. 

Computational  load  is  generally  higher  with  these  techniques  because  of  the 
eigenvalue/eigenvector  evaluations  involved.  The  particular  computer  code  used  by  the  author  is  a 
Fortran  IV  program  called  CMPLXE1G  developed  originally  at  the  University  of  Wisconsin  Computer 
Center  [8].  Offsetting  the  increased  computation,  however,  is  the  ever-increasing  speed  and  capacity  of 
modern  CPU  machines.  In  addition,  rather  efficient  eigenvalue/eigenvector  computation  methods 
[9,10]  exist  which  may  be  pursued  if  needed,  and  the  impact  of  very  large  scale  integrated  circuits 
(VLSI)  will  undoubtedly  soon  be  felt  in  modem  signal  processing. 

4.1  Snapshot  Signal  Model 

Adaptive  response  to  a  particular  signal  environment  is  based  on  algorithm-controlled  processing 
of  the  sampled  signals  from  multiple  sensors.  In  this  section  we  address  the  sampled  signals  by 
developing  in  detail  a  snapshot  signal  model  which  forms  the  basis  for  all  computer  simulations  used  in 
this  report 

Consider  a  simple  linear  array  of  K  elements  as  shown  in  Fig.  20.  The  received  signal  samples  are 
correlated  in  both  space  and  time,  giving  rise  to  a  two-dimensional  data  problem,  but  we  convert  this  to 
spatial  domain  only  by  assuming  that  narrow-band  filtering  precedes  our  spatial  domain  processing. 
Bandwidth  can  be  handled  when  necessary  via  a  spectral  line  approach  [11]  or  tapped  delay  lines  at  each 
element  [12],  but  we  did  not  consider  such  extra  complication  essential  to  the  basic  purposes  of  this 
analysis.  Thus,  the  postulated  signal  environment  on  any  given  observation  consists  of  /  narrow-band 
plane  waves  arriving  from  distinct  directions  9,.  The  RF  phase  at  the  klh  antenna  element  as  a  result 
of  the  /th  source  would  be  the  product  ottXk,  where  Xk  is  the  location  of  the  element  phase  center  with 
respect  to  the  midpoint  of  the  array  in  wavelengths,  and  at,  is  defined  as 

—  2w  sin  0/.  (1) 

This  notation  is  deliberately  chosen  to  have  the  spatial  domain  dual  of  sampling  in  the  time  domain,  so 
that  the  reader  may  readily  relate  to  the  more  familiar  spectral  analysis  variables.  Sin  9,  is  the  dual  of  a 
sinusoid  frequency  fh  and  the  Xk  locations  are  the  dual  of  time  sampling  instants  4.  Note  that  if  our 
elements  are  equally  spaced  by  a  distance  d,  then  Xk  may  be  written, 


where  X  is  the  common  RF  wavelength.  The  ratio  d/k  becomes  the  dual  of  the  sampling  time  T  with 
the  cut-off  frequency  equal  to  the  reciprocal  [13]. 
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The  complex  amplitude  of  the  l  th  source  at  the  array  midpoint  phase  center  is  p„  such  that  we  can 
now  express  the  nth  time-sampled  signal  at  the  Ath  element  as, 

Ek(n)  -  tj*(n)  +  JJ  Pi(n)gk(B,)  exp  (jot,Xk)  (3) 


where  &($/)  is  the  element  pattern  response  in  the  direction  9h  and  i j*(»)  is  the  nth  sample  from  the 
Ath  element  independent  Gaussian  receiver  noise.  (The  receiver  noise  component  is  assumed  to  be  a 
random  process  with  respect  to  both  the  time  index  n  and  the  element  index  k.)  Equation  (3)  permits 
us  to  construct  a  convenient  column  vector  of  observed  data  in  the  form. 


£(»)  —  Vp(«i)  + 1|(«) 
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1  *2JC  1 

•  VIK  _ 

(4) 


where  V  is  a  K  x  /  matrix  containing  a  column  vector  v,  for  each  of  the  /source  directions;  i.e., 

v»  -  gk(9,)  exp  (Jo>,Xk).  (5) 
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Note  that  £q.  (4)  separates  out  the  basic  variables  of  source  direction  in  the  direction  matrix  V,  source 
baseband  signal  in  the  column  vector  pin),  and  element  receiver  channel  noise  in  the  column  vector 
nin).  The  vector  E(n)  is  defined  as  the  nth  snapshot,  i.e.,  a  simultaneous  signal  sampling  across  all 
K- array  elements  at  the  nth  time  instant.  These  snapshots  would  nominally  occur  at  the  Nyquist  sam¬ 
pling  rate  corresponding  to  our  receiver  bandwidth  [13],  so  that  a  radar-oriented  person  may  view  them 
as  range  bin  time  samplings.  However,  for  source  estimation  purposes,  they  need  not  necessarily  be 
chosen  from  contiguous  range  bins  and,  in  fact,  for  most  applications  it  would  be  highly  desirable  to 
selectively  time  gate  the  snapshots  used  for  source  estimation.  For  this  simple  analysis,  let  us  postulate 
that  the  snapshots  are  gated  at  more  or  less  arbitrary  instants  of  time. 

Over  typical  processing  intervals,  the  directions  of  arrival  will  not  change  significantly,  so  that  V  is 
a  slowly  changing  matrix.  In  contrast,  tho  signals  p,in)  will  generally  vary  rapidly  with  time,  often 
unpredictably,  such  that  we  must  work  with  their  statistical  descriptions.  It  is  assumed  that  the  signals 
are  uncorrelated  with  receiver  noise.  Proceeding  then  from  Eq.  (4),  we  can  obtain  the  covariance 
matrix  R  via  application  of  the  expected  value  operator  8’,  or  ensemble  average, 

R  -  g’[E(n)E*'(n)J  (6) 

R  -  VPV *'tN  (7) 

where  N  -  &[ («)*'],  P  -  &lp(n)p(n)  *'],  *  is  the  complex  conjugate,  and  /  is  the  transpose.  N 
is  a  simple  diagonal  matrix  consisting  of  the  receiver  channel  noise  power  levels.  The  diagonal  ele¬ 
ments  of  P  represent  the  ensemble  average  power  levels  of  the  various  signal  sources,  and  off-diagonal 
elements  can  be  nonzero  if  any  correlation  exists  between  the  sources.  Note  that  correlated  far-field 
signals  can  easily  arise  if  significant  specular  reflection  or  diffraction  multipath  is  present. 

4.2  Sample  Covariance  Matrix  (SCM)  Techniques 

Antenna  array  snapshots  are  the  fundamental  input  data  to  our  processor,  whether  they  are  real 
measurements  or  computer  simulated  via  the  snapshot  signal  model  developed  in  the  previous  section. 
They  are  generally  used  to  computer  update  one  of  the  following  vectors  or  matrices  kept  in  computer 
storage: 

•  Optimum  weight  vector.  (Example— the  Howells-Applebaum  algorithm  [14]) 

•  Optimum  weight  matrix.  (Example— the  Gram-Schmidt  algorithm  [15]) 

•  Sample  covariance  matrix  (SCM).  (To  be  discussed  in  this  section) 

•  The  Inverse  of  an  SCM.  (Example— Shapard  et  al.  in  Ref.  12) 

This  section  addresses  the  SCM  in  some  detail  and,  in  particular,  describes  four  different  techniques 
that  have  been  used  in  adaptive  array  processing: 

1.  simple  block  averaging, 

2.  time  decay  averaging, 

3.  sliding-window  averaging,  and 

4.  forward-backward  subaperture  averaging. 
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The  SCM  may  be  viewed  as  a  repository  of  information  gathered  by  the  sampling  antenna  array,  in 
which  the  individual  covariance  terms  are  estimates  or  averages  derived  from  a  finite  number  of  data 
samples.  In  accordance  with  basic  statistical  theory  [16],  the  accuracy  of  the  covariance  estimates 
improves  as  the  number  of  data  samples  increases. 

4.2.1  Simple  Block  Averaging 

This  is  an  SCM  average  taken  over  N  snapshots,  and  may  be  expressed  in  the  form, 

i-  J  Ij  *<■>*<■>  *1-  (8) 

Characteristics:  storage  is  required  for  the  accumulating  SCM,  but  not  for  the  data  snapshots;  data, 
either  before  or  after  the  block  of  AT  snapshots,  are  obviously  not  included;  the  value  of  N  chosen  must 
be  large  enough  to  obtain  satisfactory  estimates  of  the  desired  processing  output  quantities;  as  each 
block  of  N  snapshots  is  completed,  M  is  passed  on  for  further  algorithm  processing  and  a  new  SCM 
accumulation  commences;  the  method  works  well  for  signal  environments  whose  statistics  change 
slowly  and  are  satisfactorily  integrated  per  choice  of  N. 

4.2.2  Time  Decay  Averaging 

This  includes  all  previous  data  in  the  SCM  average,  but  gives  ever-decreasing  weight  to  past  data; 
i.e.,  it  applies  exponential  time  decay  or  'forgetting  *  of  past  data.  Such  decay  averaging  is  the  digital 
recursive  equivalent  of  analogue  filter  integrator  circuits.  For  example,  the  equivalent  of  a  simple 
single-pole  RC  analogue  integrator  filter  [17]  may  be  expressed  in  the  recursive  backward-difference 
form, 

jR(n-l)  +  l[E(n)E(/.)*lJ  (9) 

where  —  1)  is  the  previous  stored  SCM  and  r  is  a  selectable  integration  time  constant  Charac¬ 
teristics:  storage  is  required  for  the  previous  SCM,  but  not  for  the  data  snapshots;  all  previous  data  are 
included,  but  with  exponential  time  decay  applied  to  the  older  data;  the  value  o>  r  must  be  chosen  large 
enough  to  obtain  satisfactory  integration;  R(»  —  1)  is  continuously  available  regardless  of  the  value  of 
n;  Le.,  it  is  always  current;  the  method  works  well  for  signal  environments  whose  statistics  change 
slowly,  but  may  not  be  suitable  for  abrupt  bursts  of  signal  because  of  the  inherently  long  memory 
effect 

4.2.3  Sliding-  Window  A  veraglng 

This  is  an  SCM  average  taken  over  N  snapshots,  but  here  a  window  continually  moves  with  the 
data  in  a  'bumping"  operation.  A  block  of  computer  memory  must  be  reserved  for  storage  of  the  N 
previous  snapshots, 

f - E{n  -  N) - 1 


(10) 

- E(n  -  2) - 

- E(n  -  1) - 
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When  the  current  nth  snapshot  comes  in,  the  entire  data  stack  is  bumped  up  by  1  such  that  the  top 
'oldest"  snapshot  is  bumped  off  the  stack.  In  this  method,  one  simply  keeps  the  snapshot  data  stack 
current  and  computes  the  SCM  only  when  needed, 

R0i)--Jt  I  lE(/)E(/)*'J.  (11) 

Characteristics:  storage  is  required  for  N  data  snapshots,  but  not  for  the  SCM;  after  bumping  off  of  the 
memory  stack,  old  data  is  completely  discarded;  the  value  of  N  must  be  chosen  large  enough  to  obtain 
satisfactory  integration;  the  SCM  is  computed  only  when  needed;  the  method  is  well  suited  to  signal 
environments  whose  statistics  change  rapidly  or  abruptly  because  of  the  finite-memory  stategy  and  con¬ 
tinuously  updated  snapshot  data  stack. 

4.2.4  Forward-Backward  Subaperture  Averaging 

This  is  an  excellent  technique  for  increasing  the  effective  averaging  of  our  SCM  when  needed, 
and  it  may  readily  be  implemented  if  the  antenna  array  elements  are  identical  and  equally  spaced;  Fig¬ 
ure  21  illustrates  this  technique.  We  form  a  reduced  dimension  subaperture  of  L  elements,  where  L 
must  be  less  than  the  total  number  of  array  elements  K.  Starting  from  the  left-hand  side,  the  subaper¬ 
ture  samples  its  first  snapshot  as  elements  1  through  L,  then  bumps  to  the  right  by  1  and  samples  its 
second  snapshot  as  elements  2  through  ( L  +  1),  then  bumps  to  the  right  by  1  and  samples  its  third 
snapshot  as  elements  3  through  ( L  +  2),  etc.  After  bumping  across  to  the  Ath  element,  we  will  have 
accumulated  (K  —  L  +  1)  subaperture  snapshots  from  one  overall  array  data  snapshot,  such  that  we 
can  increase  our  SCM  averaging  by  that  same  factor.  This  subaperture  motion  from  left  to  right  pro¬ 
duces  what  is  generally  termed  "forward  averaging."  The  technique  may  be  applied  to  any  SCM 
method.  For  example,  the  simple  block  averaging  of  Eq.  (8)  becomes 

»' -  jwr-i  +  i)  | "T (12> 

where  E(n,/)'-  [Et(n),  £,/+1(n),  ....  £/+i,_t(n)I  and  R/  is  the  new  reduced  L  x  L  dimension  SCM. 
Note  that  E  («,/)  may  be  expressed  as  the  matrix  product, 

E («,/)  -  I,E(n)  (13) 


where  I,  is  a  special  L  x  K  rectangular  sampling  matrix  in  which  the  /  index  denotes  the  first  column 
where  the  identity  matrix  I  begins.  For  example,  the  I,  matrix  for  L  -  3,  K  -  6,  and  i  -  2  would  be, 


I 


i 


Ij- 


-  [Oil  1 01 

010000 
001000 
0  0  0  1  0  OJ 


(14) 


An  I,  matrix  may  be  used  to  reduce  the  number  of  computations  by  multiplying  the  SCM  of  Eq.  (8)  to 
obtain 


R/“  (K-L  +  1) 


UC-L+l) 


,5 


(I/RI/J. 


(15) 


Equations  (15)  and  (12)  give  identical  results,  and  both  are  mathematical  expressions  of  the  additional 
spatial  averaging  or  "smoothing"  that  is  obtained  via  the  moving  subaperture  technique. 


Furthermore,  the  averaging  can  be  doubled  again  by  reversing  our  subaperture  at  the  right-hand 
side  and  bumping  across  to  the  left-hand  side  in  similar  fashion,  but  it  requires  conjugating  the  subaper¬ 
ture  snapshots.  This  subaperture  motion  from  right  to  left  produces  what  is  generally  termed 
"backward-averaging."  The  reader  can  verify  that  the  resulting  SCM  will  be  the  conjugate  transpose  of 
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TYPICALAPERTURE 
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BACKWARD  SHIFT 

Fit  21  —  Forward-backward  ihift  movement  for  a  reduced  dimen- 
»on  sampling  wbaperture  where  L  "  4,  along  a  linear  antenna  array 
of  K  —  15  element! 

8/.  *nd  that  we  can  combine  the  two  into  a  final  forward-backward  average  SCM  which  is  denoted  as 
the  reduced  L  x  L  matrix  Ry»: 

Ry»  -  ~  [R/  +  R/*'].  (16) 

Note  that  Ry»  is  a  symmetric  matrix,  but  is  generally  not  Toeplitz.  References  7,  18,  19  are  recom¬ 
mended  for  further  detailed  discussion  of  the  technique. 

Although  forward-backward  subaperture  averaging  is  a  very  simple  concept,  it  usually  produces 
remarkable  improvements  in  output  estimates  and,  in  addition,  becomes  crucial  to  processing  in  the  fol¬ 
lowing  situations: 

•  When  only  a  few  data  snapshots  are  available  per  SCM  computation.  Note  that  the 
method  can  be  used  even  under  the  extreme  condition  of  only  a  single  snapshot 

•  When  significant  coherence  exists  between  spatial  sources  as,  for  example,  in  multipath 
situations  involving  a  specular  reflection.  For  this  particular  condition,  the  fields  arriving 
at  the  aperture  are  nonstationary  in  space  and  the  SCM  is  not  Toeplitz  [1,5,191. 

A  caveat  concerning  this  averaging  technique  is  that,  as  the  dimension  L  of  the  subaperture 
becomes  smaller,  the  subaperture  antenna  gain,  resolution,  and  degrees  of  freedom  decrease.  Thus, 
the  advantage  of  increased  averaging  must  always  be  balanced  against  these  factors,  and  it  is  usually 
prudent  to  process  with  as  large  a  subaperture  dimension  as  possible. 

4.3  Eigenvalne/Elgenvector  Decomposition 

When  a  signal  is  known  to  consist  of  pure  sinusoids  in  white  noise,  an  appropriate  procedure  to 
find  the  unknown  frequencies  and  powers  is  the  Pisarenko  spectral-decomposition  procedure  [20], 
Although  Pisarenko’s  method  per  se  has  not  found  widespread  use,  it  has  provided  a  fundamental 
eigenanalysis  basis  for  several  closely  related  techniques  which  have  demonstrated  excellent  perfor¬ 
mance.  Among  these  are  the  algorithms  described  by  Reddi  [21],  the  MUSIC  algorithm  of  Schmidt 
[22],  the  work  of  Bienvenu  and  Kopp  [23],  the  singular  value  decomposition  (or  principal  eigenvector) 
methods  of  Kumaresan  and  Tufts  [24, 25], the  eigenassisted  method  of  Evans  et  ai.  [19],  and  the  alge¬ 
braic  approach  of  Bronez  and  Cadzow  [26].  | 

A  key  principle  in  all  of  these  techniques  is  the  geometric  vector  space  relationships  between  the 
spatial  source  vectors  and  the  eigenvectors  of  the  sample  covariance  matrix;  so  we  begin  our\  discussion 
on  this  point  From  the  theory  of  matrices,  we  know  that  a  positive  definite  Hermitian  matrix  such  as 
R  of  Eq.  (7)  can  be  diagonalized  by  a  nonsingular  orthonormal  modal  matrix  transformation  v^hich  shall 
be  defined  as  the  matrix  Q.  Furthermore,  we  know  that  the  resulting  diagonal  components  are  the 
eigenvalues  of  matrix  R.  In  accordance  with  the  usual  eigenvalue  problem  statements,  j 
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|R  ~  /3(2l|  —  0  and  Re,  —  fife,,  (17) 

the  fil  are  the  eigenvalues  (real  positive  numbers)  of  R,  I  is  the  identity  matrix,  and  e,  are  the  associ¬ 
ated  eigenvectors.  These  eigenvectors,  which  are  normalized  to  unit  Hermitian  length  and  are  orthogo¬ 
nal  to  one  another,  make  up  the  columns  of  the  Q  matrix, 


Q- 


I  t  I 

*1  *2  *3 

III 
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*K 

I 


where  e,  — 


(18) 


Diagonalization  of  R  by  the  Q  matrix  transformation  per  Eq.  (17)  may  be  written. 


[Q*'RQ]  -  \fi\)  - 


where  8tf  is  the  Kronecker  delta  symbol.  One  can  readily  show  a  construction  of  R  from  its  orthonor¬ 
mal  components,  j 

R  -  j  (20) 

Next,  we  introduce  the  term  "principal  eigenvector"  (PE)  to  mean  those  eigenvectors  which 
correspond  to  the  unique  eigenvalues  generated  by  the  spatial  source  distribution;  and  the  term  "noise 
eigenvector"  to  mean  those  eigenvectors  which  correspond  to  the  small  noise  eigenvalues  generated  by 
the  receiver  channel  noise  in  Eq.  (7).  Under  ideal  conditions,  the  noise  eigenvalues  are  all  identical 
and  equal  to  receiver  channel  noise  power  level  fil,  such  that  we  can  factor  Eq.  (20)  to  emphasize  the 
PE, 

R- ij<tf-/3.2)e,e,"‘  +  02l  (21) 

where  ;  is  the  number  of  PE.  Comparing  Eq.  (21)  with  Eq.  (7)  we  note  that  the  noise  diagonal 
matrices  are  equal;  i.e., 

fill  -  N  (22) 


fil  0  0 

0  fit  o 
0  0  fil 


fii 


(19) 
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such  that  one  may  equate  the  source  direction  vectors  with  the  PE, 

VPV'-  5}  IaI V,*'"  2J  03, 2  -  02)e,e,*'  (23) 

where  the  Ip,  I2  represents  the  expected  power  levels  of  uncorrelated  sources.  Equation  (23)  embodies 
the  key  principle  that  the  PE  are  linear  combinations  of  the  source  direction  vectors,  and  vice-versa.  In 
geometrical  language,  the  v,  define  an  /  dimensional  source  vector  space,  and  the  principal  e,  span  that 
same  vector  space.  Furthermore,  since  the  noise  eigenvectors  are  always  orthogonal  to  the  PE,  then  it 
follows  that  the  noise  eigenvectors  must  occupy  a  subspace  which  if  orthogonal  to  the  source  vector 
space.  To  put  it  another  way,  if  the  noise  eigenvectors  are  viewed  as  antenna  array  element  weights, 
then  they  should  have  pattern  nulls  at  source  direction  angles  because  of  their  orthogonality.  Despite 
the  fact  that  Eq.  (23)  is  based  on  ideal  assumptions,  it  turns  out  to  be  a  valuable  concept  for  formulat¬ 
ing  algorithms,  perhaps  because  it  is  inherently  a  noise-subtracted  relationship,  and  the  estimates  of  the 
PE  are  rather  robust 

When  working  with  finite  sets  of  data  snapshots  which  are  not  ideal,  a  nontrivial  problem  area 
arises  in  determining  which  eigenvectors  to  designate  as  principal,  and  which  ones  result  from  noise. 
This  important  problem  will  be  addressed  after  discussing  the  associated  algorithms. 

4.4  Eigenanalysis  of  Three  Algorithms 

Eigenvalue/eigenvector  decomposition  is  now  applied  to  three  different  spatial  spectrum  e.  tima- 
tion  algorithms  to  develop  the  particular  algorithm  which  was  chosen  for  processing  sample  covariance 
matrix  data.  Resultant  source  information  derived  from  these  algorithms  directly  determines  the  adap¬ 
tive  spatial  filtering  applied  prior  to  target  tracking.  These  algorithms  are 

•  the  MLM  (Maximum  Likelihood  Method), 

•  the  MUSIC  (Multiple  Signal  Classification),  and 

•  the  PEGS  (Principal  Eigenvector  Gram-Scluuidt). 

Since  these  algorithms  are  dealt  with  in  a  very  abbreviated  manner  in  this  report,  the  reader  is 
encouraged  to  consult  the  references  given  for  a  better  description  and  understanding  of  the  '  .  miques 
involved. 

4,4.1  The  MLM  (Maximum  Likelihood  Method) 

The  maximum  likelihood  spectral  estimate  is  defined  as  a  filter  designed  to  pass  the  power  in  a 
narrow  band  about  the  signal  frequency  of  interest  and  to  minimize  or  reject  all  other  frequency  com¬ 
ponents  in  an  optimal  manner  [4,27,28].  This  is  identical  to  the  use  of  a  zero-order  mainbeam  direc¬ 
tional  gain  constraint  in  adaptive  arrays  [29,30],  where  the  spatial  spectrum  would  be  estimated  by  the 
output  residual  power  P0  from  the  optimized  adapted  array  weights, 

P0(o)  -  W0*'RW0.  (24) 

Where  W0  is  the  optimum  adaptive  Wiener  filter  weight,  and 

W„-MR-,S*  (25) 

where  S*  is  the  usual  mainbeam  weight  vector  for  steering  angle  0,  and  n  is  a  complex  number.  Under 
the  zero-order  gain  constraint,  we  require  S'W„  —  1,  whereupon  ft  becomes 

ft  ■»  [S'R"IS*]-1.  (26) 
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Substituting  p.  and  W0  into  Eq.  (24)  results  in 

Po(0) 


1 

S'R-’S* 


(27) 


Upon  sweeping  the  steering  vector  S*for  a  given  covariance  matrix  inverse,  P0(9)  estimates  the  spatial 
spectrum. 


In  terms  of  eigenvalue/eigenvector  decomposition,  we  can  take  the  inverse  of  Eq.  (19)  and 
express  R-1  in  the  form, 


R"1 


•  t 


(28) 


Here,  we  see  that  this  older  algorithm  simply  uses  all  of  the  eigenvalues/eigenvectors.  One  advantage 
in  this  decomposition  is  that  Eq.  (28)  can  be  substituted  into  Eq.  (27)  to  form  a  simple  summation  of 
eigenvector  beams  referenced  to  the  receiver  noise  power  level. 


I  Po(0)  | 
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00 

0l 

g?(o) 

(29) 


where  gk(0 )  —  Sfe*.  This  permits  an  insight  into  the  peak  values  which  occur  at  the  nulls  of  the  noise 
eigenvector  beams;  i.e.,  we  get  an  evaluation  of  relative  source  power  level  if  the  source  location  is 
resolved. 


4.4.2  The  MUSIC  (Multiple  Signal  Classification) 

This  algorithm  was  suggested  by  Schmidt  {22]  to  provide  asymptotically  unbiased  estimates  of  the 
number  of  signal  sources,  directions  of  arrival,  strengths  and  cross  correlations  among  the  directional 
waveforms,  polarizations,  and  strength  of  noise/interference.  His  geometrical  vector  space  description 
and  interpretation  is  clearly  presented  and  was  used  as  the  basis  for  Section  4.3  of  this  report.  Essen¬ 
tially,  this  MUSIC  algorithm  selects  and  uses  only  the  noise  eigenvectors  to  solve  for  the  directions  of 
arrival.  This  is  tantamount  to  approximating  R-1  in  Eq.  (27)  by  the  noise  eigenvectors  only,  i.e.. 


let  R-1 


K 

I, 

(•1+1 


e,e, 


•< 


(30) 


where  q  is  the  number  of  principal  eigenvectors.  The  same  indexing  would  apply  as  noise  eigenvector 
beams  in  Eq.  (29),  where  we  note  that  the  ratio  of  eigenvalues  would  now  become  unity  (or  close  to 
it). 


This  algorithm  does  indeed  produce  very  large  peaks  in  PB(0)  for  good  covariance  matrix  esti¬ 
mates,  because  of  the  aforementioned  orthogonality  of  the  noise  eigenvectors  to  the  source  vector 
space.  Its  performance  is  usually  far  superior  to  the  older  MLM  algorithm  in  resolving  closely  spaced 
source  directions.  In  addition,  Schmidt  points  out  that  once  the  directions  of  arrival  have  been  found, 
the  direction  matrix  V,  in  Eq.  (4)  and  (7)  becomes  available  and  may  be  used  to  compute  the  source 
power  matrix  P.  We  form  the  special  matrix  U, 
where 

U  {VWV]-1V«*  (31) 


and 


UVPV'U*'-  P. 
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II'  the  matrix  U  exists  or  can  be  closely  approximated,  then  we  can  apply  it  to  the  noise-subtracted 
covariance  matrix  of  Eq.  (7)  to  solve  for  P;  i.e., 

P  -  UtR  -  N]UW  (32) 


Note  that  the  diagonal  elements  of  P  represent  power  estimates  of  the  sources,  and  that  the  nonzero 
off-diagonal  elements  represent  estimates  of  the  correlations  existing  between  partially  coherent 
sources. 

The  ability  to  solve  for  the  power  estimates  is  of  great  importance  in  distinguishing  "false  alarms” 
and  in  selecting  the  sources  of  interest  We  recommend  Refs.  19,  25,  and  26  which  are  either  related 
to  or  give  a  comparative  analysis  of  the  MUSIC  algorithm. 

4.4.3  The  PEGS  (Principal  Eigenvector  Gram-Schmidt) 

Several  eigenvalue/eigenvector  decomposition  techniques  described  in  Refs.  19,  21,  25,  and  26 
are  based  on  the  principal  eigenvectors  (PE)  together  with  some  type  of  constraint  imposed  on  the 
optimum  weight  vector. 

This  subgroup  of  PE  methods  is  of  interest  in  the  current  work  because  of  their  generally  superior 
performance  characteristics.  An  intuitive  reasoning  behind  their  use  is  that  the  estimates  of  the  PE  are 
robust;  i.e.,  they  tend  to  remain  relatively  stable  from  one  data  record  to  the  next,  whereas  the  noise 
eigenvectors  tend  to  fluctuate  because  of  noise  perturbations.  In  addition,  the  PE  methods  are 
inherently  a  noise  subtraction  technique,  similar  to  noise  power  cancellation  algorithms  which  attempt 
to  remove  the  noise  bias  term  that  appears  along  the  main  diagonal  of  the  covariance  matrix.  Direct 
application  of  such  reasoning  is  illustrated  by  the  eigenassisted  autoregressive  (EAR)  technique 
described  in  Ref.  19. 


Let  us  begin  development  of  our  PEGS  algorithm  by  decomposing  the  inverse  of  the  covariance 
matrix  as  given  in  Eq.  (28),  normalized  by  receiver  noise  power, 


Pi -Pi 


Substituting  Eq.  (33)  into  Eq.  (25)  results  in  the  optimum  Wiener  weight. 
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where 

a,  —  e*/S*, 
ft'  -  ft/pi,  and 

S*  —  a  quiescent  array  weight  vector. 

In  the  limit  of  noise-free  data,  Eq.  (34)  is  suggestive  of  a  simple  Gram-Schmidt  vector  subtraction 
from  S*  in  which  we  would  form  an  optimum  weight  that  would  be  orthogonal  to  the  PE  and,  there¬ 
fore,  orthogonal  to  the  source  direction  vectors  per  Eq.  (23).  Thus,  let  us  formulate  such  a  PEGS  algo¬ 
rithm  by  defining  the  optimum  weight  W(  from  Eq.  (34)  as, 

Wt-S*-^ja,e,  (35) 

where 


W,  possesses  the  necessary  orthogonality  to  noise-free  source  direction  vectors, 

<W„  v,>  -  0  (36) 

and  may  readily  incorporate  the  option  of  unit  Hermitian  length  if  desired, 

|W,P-1.  (37) 

The  PEGS  algorithm  as  applied  in  this  report  used  an  end-element  weighting  for  S*  i.e., 

S*-[0  0,  ....  0.  IK  (38) 


4.5  Colling  Principal  Elgent  aloes 

The  number  of  principal  eigenvalues  is  usually  directly  related  to  the  number  of  sources  which,  in 
practice,  are  not  known  and  must  be  estimated.  One  of  the  early  estimation  techniques  which  has  often 
been  used  is  the  AIC  (Akaike  information  criterion)  [31,32].  This  criterion  has  been  successfully 
applied  to  many  model  identification  problems  in  engineering  and  statistics,  including  the  well  known 
problem  of  determining  the  order  of  an  autoregressive  (AR)  process  [32].  Recent  work  reported  by 
Wax  and  Kailath  [33]  presents  a  new  approach,  based  on  the  AIC,  which  eliminates  the  need  for  any 
subjective  judgment  in  the  decision  process;  i.e.,  the  procedure  does  not  require  any  subjectively  chosen 
threshold.  This  new  approach  was  implemented  during  the  current  investigation  and  was  found  to  be 
very  effective  for  most  of  thJ  examples  tested.  In  addition  to  the  Wax-Kailath  AIC  approach,  we  also 
used  a  second  effective  technique  which  is  based  on  the  following  three  processing  operations: 

•  an  initial  decreasing-magnitude  sort, 

•  culling  per  coarse  magnitude  threshold,  and 

•  culling  per  sensitive  threshold  based  on  a  quadratic  curvefit  predictor. 

To  explain  these  operations  effectively,  a  typical  example  is  analyzed. 

Case  CSLIJIT^G  consisted  of  an  8-element  linear  array  with  elements  spaced  half-wavelength, 
receiving  signal  on  each  snapshot  from  a  10  dB  SNR  source  located  at  0s  and  signal  pulsed  on  every 
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64th  snapshot  from  a  7  dB  SNR  source  located  at  6°,  480  snapshots  of  data  used  to  form  the  sample 
covariance  matrix,  and  a  receiver  noise  power  level  of  17.64.  Thus,  the  initial  constants  were 

K  —  8  (number  of  array  elements), 

N  ■»  480  (number  of  snapshots),  and 
0 o  mm  17.64  (receiver  noise  power  level). 

Computation  of  the  eigenvalues  for  a  simple  block  average  SCM  (Section  4.2.1)  resulted  in  eight  values 
which  appeared  in  the  initial  sorting  sequence  as  1442.2,  23.8,  20.2,  19.9,  18.0,  17.5,  16.4,  and  15.7. 


The  coarse  magnitude  threshold  C0  is  computed  from  the  expression, 


C0 


(1  +  a)  0O2  + 


(39) 


where  D0  is  the  average  noise  eigenvalue  slope  and  a  is  a  coarse  magnitude  threshold  factor.  The  aver¬ 
age  noise  eigenvalue  slope  can  be  readily  determined  from  noise  eigenvalue  statistics  as  a  function  of 
the  number  of  snapshots  averaged  in  the  sample  covariance  matrix.  A  satisfactory  approximation  for 
Eq.  (39)  has  been  found  to  be 


D0 


IL 

V77 ' 


(40) 


which  results  in  D0  —  0.8  for  our  example.  The  value  of  a  must  be  chosen  large  enough  to  cover  nor¬ 
mal  receiver  noise  statistical  behavior,  and  0.24  was  found  to  be  adequate.  Substituting  these  values  for 
our  example  results  in  a  coarse  magnitude  threshold  of  C0  ”  25.9,  which  culls  out  one  principal  eigen¬ 
value  from  the  initial  sort.  Note  that  02  agrees  closely  with  the  approximation  for  the  two  closely 
spaced  sources. 


>2  ~ 


\h\2  + 


\h\2 


64 


K  +  1 


ft„. 


(41) 


Figure  22  shows  the  remaining  eigenvalues  and  the  coarse  magnitude  threshold  value.  The  reader 
can  verify  that  the  average  noise  eigenvalue  slope  is  very  close  to  the  value  given  in  Eq.  (40).  The 
number  of  eigenvalues  greater  than  C0  defines  the  initial  value  of  our  unique  eigenvalue  counter  Nu  so 
that  we  start  with  the  initial  count  of  Nu  —  1 . 


Final  culling  is  based  on  a  sensitive  magnitude  threshold  which  is  computed  from  a  quadratic 
curvefit  predictor  applied  to  the  next  largest  eigenvalue.  In  Fig.  22,  mat  next  value  is  23.8  (the  second 
eigenvalue),  and  we  use  the  remaining  six  eigenvalues  to  evaluate  the  constants  in  a  basic  quadratic 

curvefit  equation, 

/Or)  —  ax2  +  bx  +  c  (42) 


where  x  has  an  integer  range  of  —5  to  4-5.  A  sliding  window  average  is  applied  (when  sufficient  data 
values  are  available)  and  the  data  are  weighted  in  accordance  with  the  taper, 


Hx) 


1 


(43) 


based  on  an  11-point  window.  The  entire  procedure  is  listed  in  the  Appendix  as  a  Fortran  IV  computer 
code. 


A  j* 

Of  02 


Solution  of  Eq.  (42)  results  in  the  fitted  curve  shown  in  Fig.  22,  and  gives  us  the  predicted  value 
™  21.5  in  our  example.  The  sensitive  magnitude  threshold  consists  of  the  inequality, 
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Pit-  22  —  Plot  of  eigenvalues  vs  eigenvalue 
Index  for  Case  C8LIJIT4O 

This  threshold  is  exceeded  for  the  second  eigenvalue,  such  that  our  unique  eigenvalue  counter  is 
advanced  to  ’N„  —  2.  Thus,  we  determine  that  there  are  two  principal  eigenvectors  for  our  example 
case  C8L1J1T0G. 

4.6  Rootfinding  on  the  Optimum  Weights 

The  spatial  direction  location  angles  of  sources  may  be  estimated  by  applying  the  conventional 
Fourier  transform  to  evaluate  zero  locations  of  the  optimum  estimation  weights,  ns  illustrated  in  Fig.  4. 
However,  a  Fourier  transform  is  restricted  to  the  Z-plane  unit  circle  by  definition  [13],  and  may  not 
always  provide  adequate  resolution  of  roots  which  migrate  off  the  unit  circle.  A  preferable  method  is 
ihas  of  solving  for  the  exact  roots  of  the  array  polynomial.  For  our  purposes,  th  s  involves  application 
of  a  footfinding  algorithm  to  the  optimum  estimation  weights. 

- — ^  listing  of  the  Fortran  IV  computer  code  for  the  rootfinding  algorithm  used  in  this  investigation 

is  giver,  in  the  Appendix.  This  algorithm  may  be  used  to  find  the  roots  Z,  of  polynomials  of  the  form 

+  w2Z  +  w3 Z2  +  . . .  +  wKZ*~'  -  0  (45) 

where  the  wk  are  complex  array  weight  coefficients.  The  roots  are  found  by  expressing  the  polynomials 
in  terms  of  Sifjak  functions  and  by  using  the  method  of  steepest  descent  to  determine  the  zeros  [34]. 
Once  a  root  is  found,  the  polynomial  is  reduced  by  synthetic  division  and  the  process  is  repeated.  The 
algorithm  is  very  fast  and  robust. 

FicuK  Wa  illustrates  an  example  of  its  application  to  Case  C8L1J1T0C,  which  consists  of  an  8- 
element  Itaiir  a-*ray  with  elements  spaced  half-wavelength,  a  10  dB  source  located  at  -8#,  and  a  weak 
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-11  dB  source  located  at  +6°.  Simple  block  averaging  was  used  to  accumulate  a  sample  covariance 
matrix  every  480  snapshots,  and  the  PEGS  algorithm  was  then  applied  to  compute  the  optimum 
weights.  The  overlaid  plots  for  eight  successive  trials  demonstrate  the  following  points: 

•  The  strong  10  dB  source  is  located  accurately  with  very  little  variance. 

•  The  weak  -11  dB  source  shows  appreciable  variance  in  its  location  estimates,  but  note 
that  its  associated  roots  remain  close  to  the  unit  circle. 

•  All  remaining  roots  are  well  removed  from  the  unit  circle  and  may  be  disregarded,  thus 
indicating  that  only  two  sources  are  present. 

The  deviation  of  the  radial  magnitude  of  roots  from  the  Z-pIane  unit  circle  may  be  used  as  a 
coarse  culling  test  to  determine  the  number  of  sources  present.  A  threshold  value  of  0.1 5  was  selected 
for  use  with  the  PEGS  optimum  weight  roots;  i.e., 

||Z,|-  l|<  0.15.  (46) 

After  culling  the  roots  in  accordance  with  this  coarse  threshold,  the  associated  source  location  angles 
are  determined  and  used  to  construct  the  direction  matrix  V  (Section  4.1).  The  ith  root  Z,  has  the 
location  angle  9, 

Z,  -  exp  (Jn  sin  9,)  (47) 

and 

▼»  -  &(•<)  e*P  (Jo»,xk). 

where 

a*/  —  2»r  sin  9t. 

Knowing  V,  we  can  then  solve  for  the  source  power  matrix  P  as  discussed  in  Section  4.4.2,  Eqs. 
(31,32).  Figure  23(b)  illustrates  a  digital  plot  of  resultant  estimated  relative  power  with  respect  to 
receiver  noise  power  level  vs  estimated  source  location  angles  for  the  eight  successive  trials.  A  power 
level  of  -20  dB  was  selected  for  screening  out  *false  alarms"  which  might  get  through  the  root  culling, 
so  this  figure  represents  the  end  result  of  a  triple  source-screening  processing,  and  gives  us  the  final 
data: 
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4.7  Spatial  Filter  Weights 

The  target  tracking  concept  shown  in  Fig.  2  permits  great  flexibility  in  updating  the  digital  spatial 
filter  weights  because  they  are  not  developed  simply  as  a  closed-loop  feedback  adaptive  response  func¬ 
tion.  They  are  developed  primarily  on  the  basis  of  a  model  of  the  interference  environment  which  is 
updated  periodically  via  the  estimated  locations  and  relative  power  levels  of  sources  which  are  logically 
determined  to  be  interference.  In  addition,  they  can  be  modified  by  other  inputs;  for  example, 

•  the  addition  of  synthetic  sources  to  reduce  distortion  ambiguities  as  discussed  in  Section 
3.2,  or  to  counter  certain  types  of  time-modulated  interference; 


•  the  option  for  "softened*  spatial  filtering  as  a  trade-off  to  buy  less  tracking  beam  distor¬ 
tion,  as  discussed  in  this  section  and  also  Section  4.8; 

•  the  option  for  time-modulated  filter  weights,  including  range-dependent  control,  to  permit 
response  to  time/range  dependent  situations;  and 

•  the  option  for  real-time  "clean-up"  adaptive  filter  weight  perturbations.  These  are  based 
on  feedback  from  the  filter  output  residue  signals,  to  permit  fast  response  to  interference 
source  situations  which  are  changing  too  rapidly  for  the  relatively  long  periodic  updates 
from  the  spectral  estimation  algorithms. 

Therefore,  let  us  begin  by  considering  the  construction  of  a  full-aperture  dimension  estimated 
covariance  matrix  based  on  the  signal  source  data  applied  to  the  model  used  for  Eq.  (7)  in  Section  4.1 
above.  Furthermore,  let  us  solve  for  its  eigenvalues/eigenvectors  so  that  we  may  examine  the  formal 
eigenvalue/eigenvector  spatial  filter  shown  in  Fig.  24.  Here,  we  use  a  Q  matrix  beamformer  to  form 
the  set  of  orthonormal  eigenvector  beams, 

t  -  Q'E.  (48) 

Next,  let  us  introduce  an  attenuation  Ak  for  each  beam  output  such  that  the  attenuation  is  normalized 
to  receiver  noise  power  &l\  i.e.. 


where  /9a  is  the  square  root  of  £?,  and  /9*  is  the  square  root  of  the  kth  eigenvalue.  We  use  a  as  an 
attenuation  control  parameter.  When  a  —  1,  the  output  power  of  each  attenuated  beam  will  be  equal 
to  /9 because  the  eigenvector  beam  output  power  is  equal  to  its  associated  eigenvalue.  Upon  using  a 
reverse  beamformer  to  transform  back  to  element  space,  we  obtain  the  spatial  filtered  output  residue 
signal  vector, 

Ef  —  Q*ld*8u]E 

-  Q*L4*«  JQ'E.  (SO) 

Conventional  beam  weighting  S*  can  then  be  applied  to  the  filtered  output  residue  to  obtain  the  beam 
output, 

r„  -  S*'E^.  (51) 

From  Eq.  (28)  for  the  inverse  of  the  covariance  matrix,  one  can  form  an  inverse  which  is  normal¬ 
ized  to  receiver  noise  power, 

M->  -  feftjR-1  -  M  Q*',  (52) 

M-,-QL4*8wjQ*'ifa-2.  (53) 

Substituting  the  conjugate  of  Eq.  (53)  into  Eq.  (50),  we  may  rewrite  Eq.  (51)  in  the  form, 

r,-S*'M-uE  (54) 

and  then  use  Eq.  (25)  to  further  express  the  filtered  output  residue  in  terms  of  the  optimum  adaptive 
Wiener  filter  weight, 

Y0  -  W'E.  (55) 

Thus,  we  have  identical  output  from  the  optimum  adaptive  weights  and  the  formal  spatial  filter  of 
Fig.  24  when  a  —  2. 
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Fig.  24  —  Formal  eigenvalue/eigenvector 
spatial  filter  block  diagram 
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The  above  discussion  suggests  a  value  range  for  the  attenuation  control  parameter  of 

1  <  a  <  2.  (56) 

It  must  be  emphasized  that  values  of  a  less  than  2  produce  nonoptimal  filtering.  This  condition  is 
referred  to  herein  as  "softened"  spatial  filtering,  because  it  reduces  the  spatial  insertion  loss  below  the 
value  that  would  be  obtained  from  optimum  adaptive  Wiener  filter  weights.  Such  softened  insertion 
loss  filtering  is  used  as  a  trade-off  option  in  the  processing  to  buy  less  distortion  of  the  target  signals. 
Resultant  degradation  in  the  overall  SNR  is  generally  acceptable  for  most  tracking  situations,  and 
readily  reverted  to  optimal  if  needed. 

The  above  discussion  essentially  presents  two  equivalent  approaches  for  obtaining  the  spatial  filter 
weights  from  the  constructed  covariance  matrix: 

•  an  eigenvalue/eigenvector  solution  and 

•  a  matrix  inverse  solution 


Note  that  the  constructed  covariance  matrix  permits  options  such  as  adding  synthetic  sources  or  reduc¬ 
ing  or  increasing  source  power  levels.  Furthermore,  since  it  is  Toeplitz,  solutions  may  be  simplified  in 
accordance  with  readily  implemented  schemes  available  in  the  literature.  Other  computational 
economies  may  be  brought  to  bear  via  direct  use  of  the  known  array  polynomial  roots. 

4.8  Tracking  Beams  and  Distortion 


The  tracking  beams  used  throughout  this  report  are  based  on  selecting  an  adjacent  pair  of  orthog¬ 
onal,  uniform-illumination  beams  such  as  those  generated  by  a  Butler  matrix  beamformer  transforma¬ 
tion,  illustrated  in  Fig.  25.  For  example,  the  transformation  matrix  B  for  a  linear  array  with  an  even 
number  of  equally  spaced  elements,  will  have  individual  matrix  components  of  the  form 


~7W 


exp  \oi,Xk  + 


1  - 


2m 

K 


(2k-  K  -  1) 


(57) 
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where 


m  —  beam  index  (column), 
k  —  element  index  (row), 

Xk  —  element  aperture  location  from  midpoint, 
m,  -  2w  sin  9,, 

K  —  number  of  elements,  and 
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Fif.  25  —  Block  diagram  of  tracking  beam  formation 
via  Butler  matrix  beamformer 

Thus,  the  beamformer  output  vector  £  may  be  expressed  as 

£  -  B'E  (59) 

where  the  columns  of  B  are  the  individual  beamsteering  vector  weights  associated  with  each  output 
beam.  We  select  beam  number  (A72)  as  our  left  beam  of  the  pair,  and  angle  9,  as  its  pointing  direc¬ 
tion.  Note  that  the  selection  automatically  determines  the  pointing  directions  of  all  other  beams.  Fig¬ 
ure  26(a)  illustrates  the  sin  x/x  patterns  of  the  left  and  right  tracking  beams  (No.  4  and  5)  for  an  8- 
element  linear  array  with  half-wavelength  spacing.  The  pointing  direction  relationship  for  these  two 
beams  may  be  written  in  terms  of  the  Z-plane  unit  circle  circumference  segment, 

A 

1  ~it  ~  w  k*0  °(M)  ~  sin  9'^’  (6°) 

land  from  this  we  can  solve  for  the  pointing  direction  of  the  right  beam,  the  median  boresight  pointing 
direction  for  the  pair,  and  the  median  beamwidth. 
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(a)  Adjacent  pair  of  Butler  matrix  beams 
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The  familiar  sum  £  and  difference  A  tracking  beam  outputs  are  then  obtained  from  this  adjacent 
pair  via  the  3  dB  hybrid  junction  shown  in  Fig.  25. 

£- j[£m+I  +  Ej 

A- (61) 

where 

m  -  (K/2), 

and  where  we  note  that  the  difference  beam  is  in  quadrature  phase  relationship  to  the  sum  beam. 
Expressed  in  terms  of  equivalent  element  weights,  the  sum  beam  weight  S  and  the  difference  beam 
weight  D  may  be  written  via  Eqs.  (58,59) 

S  -  y  (bm+i  +  b  J 

D-y  tb„+i  -  bj.  (62) 


The  uniform  illumination  vectors  bm  and  bm+]  result  in  cosine  illuminations  for  S  and  D,  which  are 
readily  computed  for  the  b tarns  of  Fig.  26(a)  as, 


0.098 

-0.490 

0.278 

-0.416 

0.416 

-0.278 

0.490 

-0.098 

0.490 

and  D  — 

0.C98 

0.416 

0.278 

0.278 

0.416 

0.098 

0.490 

These  illuminations  are  shown  plotted  in  Fig.  26(b)  for  our  8-element  linear  array  example. 


Monopulse  tracking  [35]  involves  an  angle  estimate  for  each  pulse  (snapshot)  containing  the  tar¬ 
get,  and  it  is  computed  from  the  ratio  of  A/£.  From  the  cosine  illumination  beams  of  Fig.  26,  note 
that  we  can  form  the  approximation, 


8  is  the  angle  of  the  target  from  boresight, 

B  is  the  angle  to  the  first  null  of  the  sum  pattern,  and 
C  is  a  constant  which  depends  upon  the  particular  illuminations. 


(64) 


For  our  example,  C  ~  3.8°  and  B  “  22.0®.  Thus,  given  values  of  A /£  from  Eq.  (61),  we  can  solve  for 
the  track  angle  estimates,  8 IB,  from  Eq.  (64). 
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Next,  let  us  address  track  estimate  distortion.  Whenever  one  performs  spatial  filtering  as 
described  in  Section  4.7,  a  distortion  of  received  plane  wavefronts  occurs  because  the  spatial  insertion 
loss  generally  varies  as  a  function  of  sin  9.  This  is  shown  in  Eq.  (50)  which  relates  the  filter  output 
vector  Ey  to  the  filter  input  vector  E.  The  only  condition  under  which  E y  can  retain  the  same  phase 
slope  as  the  input  is  when  all  of  the  Ak  are  equal.  Such  distortion  can  cause  serious  tracking  errors 
from  conventional  sum  and  difference  tracking  beams. 

This  problem  was  first  addressed  in  the  literature  by  Davis,  Brennan,  and  Reed  [3],  who  proposed 
an  algorithm  for  estimating  the  angle  of  arrival,  based  on  the  outputs  of  adaptively  distorted  sum  and 
difference  beams.  They  used  approximations  to  the  optimum  angle  estimator  which  permitted  correc¬ 
tion  of  distortion  at  the  tracking  beam  boresight  position,  and  they  demonstrated  good  performance  via 
simulation  for  sidelobe  and/or  mainbeam  interference. 

The  approach  we  take  here  is  that  we  know  our  adaptive  filter  weights,  so  that  we  can  compute 
the  resultant  distortion  error  throughout  the  tracking  beam  region  and  use  it  directly  for  correction.  In 
Section  4.7  above,  we  showed  that  it  makes  no  difference  whether  we  apply  our  quiescent  beam  weights 
to  the  spatially  filtered  signals,  or  the  equivalent  adaptive  weights  to  the  unfiltered  signals.  Thus,  we 
may  apply  the  monopulse  sum  and  difference  weights  of  Eqs.  (62,  63)  to  the  spatially  filtered  output 
residue  signal  vector  of  Eqs.  (50,  51,  54,  55)  and  obtain  the  equivalent  forms. 


■S'Ey-E'M-'S-E'W, 

(65) 

D'Ey  -  E'M-‘D  -  E'Wrf, 

(66) 

where  W,  and  Wrf  are  now  the  equivalent  adapted  (and  distorted)  sum  and  difference  beam  weights. 
The  distorted  ratio  A/E  can  be  computed  for  any  direction  vector  E,  thus  giving  us  the  necessary  dis¬ 
tortion  correction  curve  across  the  entire  tracking  angle  region. 

If  there  are  no  interference  sources,  then  spatial  insertion  loss  is  zero,  there  is  no  distortion,  and 
we  have  the  ideal  linear  track  angle  estimate  case  summarized  in  Fig.  27  for  our  8-element  array.  Note 
that  Fig.  27(d)  is  a  plot  of  track  angle  estimate  8  vs  true  target  angle,  as  computed  from  A/E  values 
obtained  from  Fig.  26(b)  and  plugged  into  Eq.  (64).  Linearity  is  excellent  over  the  tracking  region  of 
±1  beamwidth.  The  abscissa  scale  in  beamwidths  is  sized  to  correspond  with  the  angle  scale  in 
degrees;  i.e.,  the  3  dB  beamwidth  of  our  8-element  array  is  about  15s. 

5.  CONCLUSIONS 

Two  closely  related  new  technologies,  spectral  estimation  techniques  and  adaptive  array  processing 
techniques,  have  been  applied  to  the  radar  systems  problem  area  of  tracking  targets  in  the  presence  of 
nearby  interference  sources.  It  appears  that  the  combination  of  these  technologies  does  indeed  offer 
some  attractive  possible  solutions.  This  application  area  was  addressed  via  an  all-digital  receive  system 
concept  whereby  adaptive  spatial  filter  weights  are  periodically  updated  on  the  basis  of  source  estimates, 
and  the  filter  output  residue  signals  are  then  sifted  for  targets.  Several  simulation  examples  were 
included  to  illustrate  some  of  the  processing  step  outputs,  "softened”  spatial  filter  insertion  loss,  track¬ 
ing  beam  distortions,  correction  for  wavefront  distortions,  and  satisfactory  tracking  in  the  presence  of 
strong  multiple  interference  sources  located  within  a  beamwidth  of  the  target.  Future  plans  include 
further  evaluation  of  performance  characteristics  via  mere  stringent/realistic  computer  simulations  plus 
testing  with  experimental  data. 
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Appendix 

FORTRAN  IV  COMPUTER  CODE  LISTINGS 


1.  DETERMINATION  OF  NUMBER  OF  UNIQUE  EIGENVALUES 

C  DETERMINATION  OF  NUMBER  OF  UNIQUE  EIGENVALUES 
C  FIRST  APPLV  COARSE  MAGNITUDE  LIMIT 
10*0  (FLAG  FOR  L0 
L»*e  | NOISE  EIGENVALUE  COUNTER 
DO  74  N-l.NEO 
IF(£IG(N).GT.BL0)GO  TO  74 
L0*L0+1 
74  CONTINUE 

IF(HOD<IPO,2>.EO.O)GO  TO  88 
76  IF(L0.LT.2)GO  TO  88  (SKIP  DIFF.  LOOPS 
C  NEXT  DETERMINE  DEVIATION  FROM  CURVEFIT 

J0-L0-U  (NUMBER  OF  CURVEFITS  ALLOWED 
IF( J0.GT.0)GO  TO  78 
JO*  1 

78  DO  86  J-1.J0 
JRITE(4)L0 
S0-0.0 
DR2-0.0 
DO  80  L— S.4 
N-L0-6+L+J 
IF CN.LT.3 )GO  TO  80 
M-NE0+1-N 

X0*EIG(H)-2.0*EIG<HM  >+EIG<R*2> 

DR2*DR2+X0*TA0( L ) 

S0*S0+TA0(L) 

80  CONTINUE 

DR2-DR2/S0  , AVERAGED  VALUE  OF  2ND  DIFF. 

IF (DR2.LT . -*(0.32<DF0 )  )DR2*-(0. 32JDF0  )  , LIMIT  NEC.  DR2 

IF(DR2.GT.+(0.16*DF0)iDR2*+<e.l6SrF0)  , LIMIT  POS.  DR2 

DR1-0.0 
S0-0.0 

DO  82  L—5.S  (SLIDING  UINDOU  SAMPLING 

N-L0-6+L+J 

IF(N.LT.2)G0  TO  82 

M-NE0*1-N 

X0-EIGCH )-EIG(H+l ) 

XFIX0.GT. O.OXDF0 ) )X0*  O.OtDFO )  (THROW  OUT  LARGE  DIFF. 
DRl*DRl+(X0-DR3tFL0AT(L  > )*TA0<L ) 

S0*S0+TA0(L )  (SUM  OF  TAPERS 
82  CONTINUE 

D81*DRt'S8  (AVERAGED  VALUE  OF  1ST  DIFF.- 
C  COMPARE  PREDICTED  VALUE 
A0-O.0 
DO  SN  L—5,4 
N-L0-6+X+J 
IF(N.LT.1)G0  TO  84 
M-ME0+1-N 
XO-FLOAT(L) 

A3*A0+EIC(H  )-DRi*X0- (DR2/2.0 )»<X0**2 ) 

84  CONTINUE 

A0-A0/FLOATIL0-1  ) 

G0-A0+DRl*S.0+(DR2/'2. 0)125.0  (PREDICT  L*S 
K-NE0+J-LO 

IF( ( EIG(K  )-G0 ) .GT. (2.0SDF0 ) )I0*1  (EIGENVALUE  TEST  FLAG 
WRITE ( 4. 12  )I0.DR2,DR1. AO.GO.EIGCK) 

86  CONTINUE  (END  OF  J  LOOP 
IF ( I0.EQ.O )C0  TO  88 
L0*L0-1  (REDUCE  L0 
10*0  (RESET  10 
CO  TO  76 

88  NU0*NE0-L0-MU00  (UNIQUE  EIGENVALUES 
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2.  ROOTFINDING  ROUTINE 


c 


c 


c 


ROOT  FINDING  ROUTINE 
609  00  692  N-1.NR0 

UTE(NiaCHPLX< 999. 999.999.999) 

692  CONTINUE 

N0*NR0  , CURRENT  ROOT  INDEX 
IFCN0.EQ. 1  )G0  TO  644 

619  2-CNPLX(0. 1.1.9) 

SFXC  0  )*CMPLX(  1 .0,9.6 ) 

5FX( 1  )*2 

L0-9  j ITERATION  COUNTER  INDEX 
INSERT  ESTIMATED  ROOTS  AT  THIS  POINT 
IFOCR0.EQ.0  )C0  TO  612 
IF( IQ9 . GT .KDF  )G0  TO  C12 
A9*PI0<SIN(PI0tDF0( IQ9)/189.0 ) 
X*COS(A0) 

V*SIN(Ae) 

2*CMPLX(X,V) 

SFXC 1  )*2 
109*109+1 
612  CONTINUE 
1.0 

GO  TO  629  j INITIAL  CONFUTE  FC2 ) 
614  1*1 
G0-F0 
H»0 

ALF*CNPLX(0. 0,0.0) 

L0-L0+1  ; ADUANCE  ITERATION  COUNT 

DO  616  K*1,N0 

ALF* ALF +FLOAT ( K  )XCFE  C  K  )*SFX C K-l ) 
616  CONTINUE 

A0*CABS( ALF  )tS2 
DIF— BETJCONJGC  ALF  )/A0 
618  N-N+l 

SFX(1)*2+DIF 

620  CONTINUE 
COHPUTE  UALUE  OF  F<2) 


642  CONTINUE 

N0-N0-1  , REDUCE  N0  INDEX 

IF(N9.NE.l )G0  TO  610  jTEST  FOR  LAST  ROOT 
844  BET*CFEC 1 ) 

A9*CABS(BET)«2 

UTEC1 )— CFEC9)tCONJGCBET)/A0 

GO  TO  250  jROOTS  COMPLETED.  RETURN  TO  IBN  LOOP 


A9aCABSCSFXCl ) )*X2 
8aB.9(REALCSFXC 1 ) ) 

N8-N0-2 
DO  624  K*0,N2 

SFXCK+2)aBSSFXCK+l )-A9tSFX(K ) 

624  CONTINUE 

1ET-CNPLXC0. 6,9.0) 

DO  626  Ka0,N0 
1ET*BET+CFECKXSFXCK) 

626  CONTINUE 

F9aCABS(BET) 

IFCF9.LT.TL1 )G0  TO  640  ;FCZ>  NULL  TEST 
IFCI .E0.6 )G0  TO  614 
IFCF9.GE.G9 )G0  TO  630 

IF (CAB  CDIF ) . LT.TL2 )C0  TO  640  jROOT  TEST 

1FCL9.GT . IT9)G0  TO  638 

Z-SFXCl) 

GO  TO  614  jNEXT  ITERATION 

639  IF(H.CT.29)G0  TO  632 
DIF-DIF/4.0 

GO  TO  618 

632  LFlal 

GO  TO  640 
638  LF2*1 

640  CONTINUE 

C  ROOT  FOUND.  STORE  AND  SUBTRACT  IT  OUT 
UTECN0  )*SFX( 1 ) 

ALF*CFECN9> 

CFECN0)*CNPLX(9. 0,0.0) 

Z*SFXC1) 

N1-N0-1 

DO  642  Ja9,Nl 

K-NI-J 

COP-CFECIC)  jSAUE  CFECK ) 

CFECK  )aALF+Z*CFECIC+l ) 

ALF-COP 
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